#!/bin/bash -e

# snp for parallel using torque's depend

function info() {
echo Usage: `basename $0` '[-lnsgGace] reads1.fq.gz reads2.fq.gz'
exit 65
}

while getopts  ":l:p:f:n:s:gGa:ce:" opts; do
	case  $opts  in
	a) cut=$OPTARG;;
	c) cnv=T;;
	l) interval=$OPTARG;;
	n) nums=$OPTARG;;
	s) sample_name=$OPTARG;;
	g) single_cell=T;;
	G) somatic=T;;
	p) out_prefix=$OPTARG;;
	f) suffix=$OPTARG;;
	e) num_pre=$OPTARG;;
	*) info;;
	esac
done
shift $(($OPTIND - 1))


if [ $# -lt 2 ]; then info; fi



r1=`rel_abs.sh $1`
r2=`rel_abs.sh $2` # must run before mkdir=true

makedir=true

if test -z "$interval"; then interval=$data_path/intervals/1/b37.bed; fi
if test -z "$num_pre"; then num_pre=400; fi

. $var


split_fastq=split_fastq.sh

if test -z "$cut"; then
cutadp=cutadp.sh
elif test "$cut" = "general"; then
cutadp=cutadp.sh
elif test "$cut" = "malbac"; then
cutadp=cutadp_malbac.sh
elif test "$cut" = "ctc80"; then
cutadp=cutadp_ctc80_discard_trimmed.sh
fi

aln=aln_bwa_bt.sh
sam2b=sam2b_samtools.sh
merge=merge_bam.sh
hmmcopy=hmmcopyb38.sh
split_bam=split_bam_0701.sh
bed_split=bed-split.py
extract_bam=extract_bam_bed.sh
extract_bam_bed=extract_bam_bed.sh
markdup=deldup.sh
realn=realign_bt.sh
recal=recal0630.sh
ug=ug.sh
annos=annos0711.sh
combine_vcfs=combine_vcfs.sh
impact=impact.sh
format_simple=format_simple_0629.sh

fq_rename=snpp_fq_rename.sh


cpu2=2
cpu4=4
cpu8=8
mem4=4
mem8=8
mem16=16
mem32=32


fq_list=$out_prefix.fq.list.txt
cutadp_list=$out_prefix.cutadp.list.txt
sam_list=$out_prefix.sam.list.txt
bam_list=$out_prefix.bam.list.txt
bed_list=$out_prefix.bed.list.txt
merge_list=$out_prefix.merge.list.txt
split_list=$out_prefix.split.list.txt
markdup_list=$out_prefix.markdup.list.txt

i_bed_list=$out_prefix.i.bed.list.txt
l_bed_list=$out_prefix.l.bed.list.txt
realn_list=$out_prefix.realn.list.txt
recal_list=$out_prefix.recal.list.txt
ug_list=$out_prefix.ug.list.txt
high_list=$out_prefix.high.list.txt
low_list=$out_prefix.low.list.txt
hl_list=$out_prefix.hl.list.txt
combine_list=$out_prefix.combine.list.txt
impact_high_list=$out_prefix.impact.high.list.txt
impact_low_list=$out_prefix.impact.low.list.txt
impact_hl_list=$out_prefix.impact.hl.list.txt
end_list=$out_prefix.end.list.txt


# split1_id=$(sub -d"synccount:1" -F -nsplit -m2 -p2 "pigz -cdf $r1|split_fastq.sh -l4000000 -p$out_prefix.1.split -")

# split2_id=$(sub -d"syncwith:$split1_id" -F -nsplit -m2 -p2 "pigz -cdf $r2|split_fastq.sh -l4000000 -p$out_prefix.2.split -")

split1_id=$(sub -F -nsplt1 -m2 -p2 "pigz -cdf $r1|split_fastq.sh -l4000000 -p$out_prefix.1.split -")

split2_id=$(sub -F -nsplt2 -m2 -p2 "pigz -cdf $r2|split_fastq.sh -l4000000 -p$out_prefix.2.split -")


fq_rename_id=$(sub -d"afterok:$split1_id:$split2_id" -F -nrenam -m1 -p1 $fq_rename)


for i in $(seq $num_pre); do

cutadp_id[$i]=$(sub -d"afterok:$fq_rename_id" -F -nadp.$i -m2 -p1 $cutadp -p$out_prefix.$i $out_prefix.1.split.$i.fastq $out_prefix.2.split.$i.fastq)

aln_id[$i]=$(sub -d"afterok:${cutadp_id[$i]}" -F -naln.$i -m6 -p4 $aln -p$out_prefix.$i -s$sample_name -t4 $out_prefix.$i.1.cut.fastq $out_prefix.$i.2.cut.fastq)

sam2b_id[$i]=$(sub -d"afterok:${aln_id[$i]}" -F -n2b.$i -m4 -p2 $sam2b -p$out_prefix.$i -t2 $out_prefix.$i.sam)

sam2b_ids=$sam2b_ids:${sam2b_id[$i]}
done

merge_id=$(sub -d"afterany$sam2b_ids" -F -nmerg -m16 -p6 $merge -t6 $out_prefix.*.sort.bam)

if test -z "$nums"; then nums=120; fi

cat $interval|cut -f1-3 >$out_prefix.bed

$bed_split $out_prefix.bed $nums

for i in $(seq $nums); do

extra_id[$i]=$(sub -d"afterok:$merge_id" -F -nextra.$i -m4 -p2 $extract_bam -p$out_prefix.$i $out_prefix.merge.bam $out_prefix.$i.bed)

mark_id[$i]=$(sub -d"afterok:${extra_id[$i]}" -F -nmark.$i -m4 -p1 $markdup -p$out_prefix.$i -m $out_prefix.$i.extract.bam)

realn_id[$i]=$(sub -d"afterok:${mark_id[$i]}" -F -nrealn.$i -m8 -p4 $realn -p$out_prefix.$i $out_prefix.$i.markdup.bam)

recal_id[$i]=$(sub -d"afterok:${realn_id[$i]}" -F -nrecal.$i -m8 -p4 $recal -p$out_prefix.$i $out_prefix.$i.realn.bam)

ug_id[$i]=$(sub -d"afterok:${recal_id[$i]}" -F -nug.$i -m12 -p6 $ug -p$out_prefix.$i $out_prefix.$i.realn.recal.bam)

ano_id[$i]=$(sub -d"afterok:${ug_id[$i]}" -F -nano.$i -m12 -p6 $annos -p$out_prefix.$i $out_prefix.$i.snp.vcf $out_prefix.$i.indel.vcf )

echo $out_prefix.$i.high.vcf >> $high_list
echo $out_prefix.$i.low.vcf  >> $low_list
ano_ids=$ano_ids:${ano_id[$i]}
done

combh_id=$(sub -d"afterok$ano_ids" -F -ncombh -m4 -p2 $combine_vcfs -p$out_prefix.high $high_list)

combl_id=$(sub -d"afterok$ano_ids" -F -ncombl -m4 -p2 $combine_vcfs -p$out_prefix.low $low_list)


imph_id=$(sub -d"afterok:$combh_id" -F -nimph -m4 -p2 $impact $out_prefix.high.combine.vcf)

impl_id=$(sub -d"afterok:$combl_id" -F -nimpl -m4 -p2 $impact $out_prefix.low.combine.vcf)


fmth_id=$(sub -d"afterok:$imph_id" -F -nfmth -m10 -p1 $format_simple -p$out_prefix.high $out_prefix.high.combine.af.hg19_multianno.vcf)

fmtl_id=$(sub -d"afterok:$impl_id" -F -nfmtl -m10 -p1 $format_simple -p$out_prefix.low $out_prefix.low.combine.af.hg19_multianno.vcf)

end_id=$(sub -d"afterok:$fmth_id:$fmtl_id" -F -nen -m1 -p1 $end `basename $0` $@)

. $cmd_done